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ABSTRACT 

The arrival directions of multi-TeV cosmic rays show significant anisotropies at small angular scales. 
It has been argued that this small-scale structure can naturally arise from cosmic ray scattering in 
local turbulent magnetic fields that distort a global dipole anisotropy set by diffusion. We study this 
effect in terms of the power spectrum of cosmic ray arrival directions and show that the strength of 
small-scale anisotropies is related to properties of relative diffusion. We provide a formalism for how 
these power spectra can be inferred from simulations and motivate a simple analytic extension of the 
ensemble-averaged diffusion equation that can account for the effect. 

Subject headings: cosmic rays, diffusion, magnetic fields, turbulence 


1. INTRODUCTION 

The propagation of cosmic rays (CRs) through turbu¬ 
lent magnetic fields is effectively described as a diffusion 
process. This approach provides an excellent descrip¬ 
tion of the observed spectra and chemical abundances 
of Galactic CRs up the CR knee at a few PeV. The ar¬ 
rival direction of CRs is mostly isotropized in this process 
and the location of the sources is hidden. The first order 
correction to the isotropic distribution is a small dipole 
anisotropy that scales with the local CR gradient accord¬ 
ing to Fick’s law. The phase and strength of this dipole 
is expected to be a combined effect of the relative mo¬ 
tion of the solar system with respect to the frame where 
CRs are isotropic (Compton and Getting 1935), the den¬ 
sity gradient of CRs, the strength of diffusion, and the 
strength and orientation of regular magnetic fields (Er- 
lykin and Wolfendale 2006; Blasi and Amato 2012; Pohl 
and Eichler 2013; Mertsch and Funk 2015; Schwadron 
et al. 2014). 

Various CR observatories have revealed anisotropies 
in the arrival directions of TeV- PeV Galactic CRs on 
large and small angular scales (Amenomori et al. 2005; 
Amenomori 2006; Guillian et al. 2007; Abdo et al. 2008, 
2009; Bartoli et al. 2013; Aglietta et al. 2009; Abbasi 
et al. 2011; Aartsen et al. 2013; Abeysekara et al. 2014). 
In this Letter we will focus on small-scale anisotropies 
that have been studied via a power spectrum analysis 
by IceCube (Aartsen et al. 2013) and HAWC (Abey¬ 
sekara et al. 2014) (see Fig. 2). Giacinti and Sigl (2012) 
have argued that small-scale anisotropies can naturally 
emerge as a next-to-leading order correction of the dif¬ 
fusive approximation if one takes into account the local 
scattering of CRs within the diffusive scattering length. 
Ahlers (2014) showed that under certain simplifying as¬ 
sumptions these small-scale fluctuations necessarily arise 
from the application of Liouville’s theorem. Assuming a 
simple hierarchical evolution equation of the multipoles 
predicts a power spectrum that matches well the obser¬ 
vations. 

Fluctuations of the arrival direction of CRs induced by 


local turbulent fields are not directly accessible in stan¬ 
dard diffusion theory. This becomes obvious if one con¬ 
siders the angular power spectrum of the CR phase-space 
distribution (PSD) f(t, r, p) defined as 

Q = ^ J d Pi /" d P2-P^(PiP2)/i/2, (1) 

where p j denotes a unit vector of p, and Pg are Legen¬ 
dre polynomials of degree l. Here and in the following we 
use the abbreviation = /(0, 0, p.,) for the local PSD. 
In general, the small-scale magnetic field that determines 
the PSD is unknown. Therefore, the PSD can only be 
predicted as an average over a statistical ensemble of ran¬ 
dom small-scale magnetic fields, usually characterized by 
its power spectrum. In ergodic systems, the average over 
turbulent magnetic field configurations, be., the ensem¬ 
ble average denoted in the following by (...), is equivalent 
to a time-average. However, the ensemble average of (Cf) 
is proportional to (/ 1 / 2 ), which is larger than the prod¬ 
uct of the ensemble averages, be., the standard solution 
of diffusion theory. 

2. RELATIVE DIFFUSION 

On large timescales, the motion of CRs in turbu¬ 
lent magnetic fields can be effectively approximated as 
a diffusion process. The angular-integrated PSD n = 
n(t , r) = / dp f(t, r, p) then satisfies a diffusion equa¬ 
tion, dtn ~ V(KVn), where the diffusion tensor K can 
be expressed as 

B; Jj, Sij BiBj e ljk B k . 

A « = T^r + ^^ + v7r- (2) 

Here, B denotes a unit vector along the regular magnetic 
field. In general, the (scattering) rates mi, and va are 
independent of one another. The presence of a small local 
gradient Vn induces a dipole anisotropy in the quasi¬ 
stationary solution as a result of Fick’s law, 47r(/) ~ n — 
3pKVn. Therefore, standard diffusion theory predicts 
that the dipole strength of the relative intensity SI = 
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( 47 t/ — n)/n is given as 
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We will see in the following that this relation becomes 
modified once we consider corrections of the PSD prod¬ 
uct in the ensemble average. This will also introduce 
multipoles at small angular scales, which can be related 
to properties of relative diffusion. 

We now study the effect of small local fluctuations of 
the PSD around the ensemble average, Sf = f — (/). 
According to Liouville’s theorem we can relate the local 
(be., r = 0) PSD f, = /(0,0,p,) to the contribution 
backtracked along CR trajectories to an arbitrary time, 


47 Tfi ~ 47r5/(-r,r < (-T),p i (-T)) 

+ n+[r i (-r)-3p i (-T)K]Vn, (4) 

where n and S7n denote the local CR density and gradi¬ 
ent and r %{—T) and pj(— T) are the position and momen¬ 
tum of a CR (that is at position r* = 0 and pj(0) = p; 
at time t = 0) at backtracking time T. Now, in the 
limit of large T the last term in Eq. (4) is dominated 
by the third term scaling with the position of the par¬ 
ticle. Also, for two momenta pi ^ p -2 we can assume 
that the ensemble average of fluctuations are uncorre¬ 
lated, (5/i(—T)<5/ 2 (—T)) ~ 0, for sufficiently large back¬ 
tracking times when the CR trajectories eventually sep¬ 
arate. In the degenerate case pi = P 2 the two back¬ 
tracked CR trajectories stay correlated over arbitrarily 
long backtracking times. It will be sufficient to assume 
that T)) 2 ) remains finite. We can then express 

the multipole spectrum of the ensemble-averaged rela¬ 
tive intensity as the limit 
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(5) 


di being shorthand for d/da 
Note that the £ > 1 multipole spectrum is generated 
through relative diffusion: it can be easily seen that the 
sum over all ensemble-averaged multipoles of the relative 
intensity can be expressed via the symmetric part of the 
diffusion tensor (ri(—T)rj(—T)) —> 27TQ) in the limit of 
large backtracking times T, 


a 2 = 1, r L /L c = 0.1, WA = 0.01, KJL C = 100 



Figure 1. Evolution of the ensemble-averaged power spectrum 
(5) for a CR gradient parallel (solid lines) and perpendicular (dot¬ 
ted lines) to the regular magnetic field and the 3D turbulence 
model discussed in the main text. We show results in terms of 
the dipole (Ci) (black), monopole (Co) (blue), and medium-^ mul¬ 
tipoles (green). We also show the asymptotic noise level (9) (red) 
and the dipole prediction (3) of standard diffusion (magenta) eval¬ 
uated by the replacement (ruT 2 j) —> ( fn)(r 2 j) in Eq. (5). 

with Ari 2 = ri — r 2 - Therefore, the sum of multipoles 
£ > 1 is related to the relative diffusion tensor. For un¬ 
correlated particle trajectories, this expression reduces 
to the normal diffusion tensor. However, particle tra¬ 
jectories with a small relative opening angle will follow 
similar trajectories and the relative contribution ( 8 ) re¬ 
mains small over long timescales. Note that the multi¬ 
poles in Eq. (5) are expected to be finite in the limit of 
large backtracking times since particle trajectories with 
arbitrarily small opening angles will eventually become 
uncorrelated, (ru(—T)r 2 j{—T)) — > 0 . 


T£(2* + l)(C,)(T)^2TJ7jAA5. (6) 

^>0 

On the other hand, the average monopole contribution 
in this limit can be expressed as 

A(Co>(r) = 2r(^-xj.)^, (7) 

where the symmetric part of the relative diffusion tensor 
is defined as 


K h = 


dpi 

4-7T 
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3. SIMULATION 

In the following, we will study the development of 
small-scale anisotropies via numerical simulations (see 
also Giacinti and Sigl (2012); Ahlers (2015); Lopez- 
Barquero (2015); Rettig (2015)). We follow the approach 
of Giacalone and Jokipii (1999) and define a three- 
dimensional (3D) turbulent magnetic field as the sum 
dB(i-) = £„=i <5B„ cos(k„r +/3 n ) with N random phases 
P n and wave vectors k„ with 3D random orientations, on 
top of a regular field B 0 . The wave vector amplitudes 
k n range from k m - ln to fc max with equal logarithmic steps. 
The vectors dB„ have 3D random orientations subject to 
the conditions dB„ _L k„ and |<5B„| 2 oc kn/ (1 + ( k n L c ) 7 ) 
with a coherence scale L c . We assume a Kolmogorov- 
type phenomenology with 7 = 11/3 and the strength of 
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the turbulence relative to the regular magnetic field is 
fixed via a parameter a 2 and the normalization of the 
spatial average <5B 2 = cr 2 B 2 . 

For the Galactic environment we expect a regular mag¬ 
netic fields of the order of 6 pG and a turbulent compo¬ 
nent with a 2 < 1 and a coherence scale L c ~ 100 pc. 
The Larmor radius of 1-10 TeV CRs is then of the order 
of Tl — 10~ 5 L C corresponding to a few hundred AU. A 
numerical study of this configuration with the methods 
described above is very time-consuming since the asymp¬ 
totic behavior relevant for the multipole configuration 
is reached at a very late stage. Instead, here we will 
study a simple turbulent configuration that allows us to 
reach the necessary statistics to study the predicted ef¬ 
fect qualitatively. We choose cr 2 = 1 with A m i n = 0.01L C 
and A max = 100L C with k = 2n/X and fix the rigidity 
of the CRs to tl = 0.1 L c . We sample over 120 dif¬ 
ferent turbulent magnetic field configurations. For each 
field configuration we uniformly sample 12288 CR ori¬ 
entations p, (0) following the HEALPix parametrization 
(n s ide = 32) (Gorski et al. 2005), which we backtrack 
through the static magnetic field B(r) = B^e, + <5B(r) 
to find the initial position rj(—T). 

Figure 1 shows the power spectrum at different back¬ 
tracking times determined via a multipole expansion 
of the relative intensity sky maps using the HEALPix 
utilities. We show the cases of CR gradients parallel 
(solid lines) and perpendicular (dotted lines) to B 0 . The 
black lines show the evolution of the ensemble-averaged 
dipole induced by the CR gradient. The magenta lines 
correspond to the standard dipole anisotropy defined 
via Eq. (3) evaluated by the replacement (ruX 2 j) —► 
(rii){r 2 j) in Eq. (5). As expected, this dipole estimate is 
smaller than the ensemble average. Note that for the case 
of perpendicular diffusion (dotted lines) the difference in 
the asymptotic values is about one order of magnitude. 
The green lines show the multipole power for 2 < £ < 9. 
We also show the expected monopole component (7) as 
blue lines. At large times, the high-<? components in the 
multipole expansion are dominated by noise, due to the 
finite number of trajectories. The noise level can be esti¬ 
mated for large backtracking times T as pixel shot noise 


AT 


4?r ot rys dindjU 

n- 2TK '’^T- 


(9) 


This is indicated as red lines in Fig. 1 and clearly influ¬ 
ences the level of high-f multipoles (green lines) in the 
map. 

The best estimator for the true power spectrum is then 
Ci = (Ci) — Af and the variance (excluding cosmic vari¬ 
ance) is given by var(C^) ~ 2Af 2 /(2£+l); see,e.<?., Camp¬ 
bell (2015). In Fig. 2 we show the estimators Ci normal¬ 
ized to the dipole C\ at a backtracking time f IT = 100, 
corresponding to the timescale of the transition into the 
diffusion regime (cf. Fig. 1). Given the finite number of 
trajectories, it might be difficult to disentangle noise vari¬ 
ance and cosmic variance in the simulations; we hence 
limit ourselves to the former only. Similar to the sim¬ 
ulation, the sensitivity of experimental data to high-t' 
multipoles is limited to shot noise indicated as dotted 
lines in the plot. Overall, the simulation agrees well with 
the IceCube and HAWC data. For an isotropic Gaus¬ 


0 2 = 1, T L /L c = 0.1, Kmn/Lc = 0.01, Vax/A = 100, SIT = 100 



Figure 2. Power spectrum estimator ( A; = (Op) — J\f (normal¬ 
ized to Ci) for the parallel (filled red diamonds) and perpendic¬ 
ular (filled green triangles) CR gradient for the same simulations 
as shown in Fig. 1 at a backtracking time SIT = 100. We also 
show the IceCube (Aartsen et al. 2013) (blue open circles) and the 
HAWC (Abeysekara et al. 2014) (magenta open squares) power 
spectra renormalized by a factor of 4 X 10® for better compari¬ 
son. The different noise level of the data is indicated as dotted 
lines. The three gray lines correspond to the prediction of a rela¬ 
tive scattering term v r (x) oc (1 — x) p in Eq. (15) for three different 
values of p. 

sian field, deviations at multipoles of £ < 5 can be due 
to cosmic variance with var(Cjj) ~ 2C 2 /(2£+ 1). Fur¬ 
thermore, the limited field of view of the observatories 
introduces additional systematic errors, in particular for 
low-f multipoles, which are not included in the error bars 
of Fig. 2. 


4. GENERALIZED BGK FORMALISM 

In the following, we derive the asymptotic form of the 
power spectrum via a microscopic description of relative 
diffusion. As before, we split the PSD into the aver¬ 
age and fluctuating parts, /=(/) + Sf. The ensemble- 
averaged Boltzmann equation can be written as (Jones 
1990) 

dt(f) + pV(/) - znL(/> = (iuUf) . (10) 

Here we have introduced the angular momentum op¬ 
erator L a = —ieabcPbd/dp c and the rotation vectors 
f2 = cBq/po and u> = eSB/po- 

The rhs of Eq. (10) encodes the particle scattering in 
the random fields. The influence of the turbulence is 
typically approximated as a friction term introduced by 
Bhatnagar et al. (1954) (BGK). In its original version 
this term drives the ensemble-averaged distribution (/) 
to its isotropic mean n with an effective relaxation rate v, 
i.e., (*u>L<5/) ~ —i/((/) — n/(47r)). The parameters v\\, 
and va are the effective parallel, perpendicular, and 
axial scattering rates, and can be related to the BGK rate 
v via F|| = v, iz_l/ rx|| = 1 + fl 2 /^ 2 , and va/^ = 1 + ^y/D 2 . 

Substituting into Eq. (2) we find that the standard cf 
and Ci (3) scale as 1 /fm and 1 /iz_l, respectively, similar 
to the sum (6) in the relative diffusion approach. 

We can generalize the BGK approximation of the fric¬ 
tion term to higher multipoles if we consider the effect 
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of turbulence on the particle momenta as a Brownian 
motion on a sphere (Yosida 1949). For a distribution 
/ on a sphere the Brownian diffusion of momenta is 
described via dtf = — (v/2)'L 2 f. In particular, if ini¬ 
tially the configuration is localized, /(0, p) = (5(pp 0 — 1), 
then the solution at later times is given as f(t, p) = 

Sfm y ?m(p) 5 ^m(Po) e ^ (f+1)W/2 - This reduces to a two- 
dimensional Gaussian distribution of the opening angle 
if) <C 1 defined via cos ij) = pp 0 with a width a 2 ~ vt. 
Therefore, for the rhs of Eq. (10) we make the ansatz 

<*«L/> ~ -(v/2)L 2 (f), (11) 

which is a natural generalization of the BGK approach. 
Note that this simple treatment of the scattering term is 
only approximate. Depending on the strength and coher¬ 
ence length of the turbulent magnetic field, the effective 
relaxation can in general be anisotropic even for a fully 
isotropic turbulence. However, in our case we are in¬ 
terested in the qualitative behavior and development of 
high-£ multipoles, and the first order BGK approxima¬ 
tion seems appropriate. 

Now, for the asymptotic solution of (Cg) at large back¬ 
tracking times we have to look for stationary solutions of 
the equation 

dt(hh) = (fi (-p 2 V + iu>2L + *OL) / 2 ) + (10 2). 

( 12 ) 

As before, we make the ansatz /) = (ff) + Sfi with the 
stationary solution (fi). Assuming that the spatial de¬ 
pendence of Sfi is small compared to (fi) we can re¬ 
place the gradient term in Eq. (12) via (/ip 2 V/ 2 ) — 
—3/(47r) 2 piV?rp2KV?r. On the other hand, the second 
and third terms in Eq. (12) are now non-trivial terms 
in the evolution that will lead to mixing and damping 
of multipoles. Analogous to the BGK-type approxima¬ 
tion (Bhatnagar et al. 1954) of Eq. (11), we make the 
ansatz, 


((iu>i Li + «£a>2L 2 )/i/ 2 ) 

M)— 


L| 


fx)3 2 


(/ 1 / 2 ), (13) 


with x = P 1 P 2 and J = Li + L 2 . The relative and cor¬ 
related scattering rates, v r and v c , respectively, depend 
on the relative distance of the trajectories at early times, 
which in turn depends on the trajectories’ opening an¬ 
gle. Note that Eq. (13) is necessarily symmetric under 
the interchange of pi •<-> p 2 as well as p, O p; and 
is the most general linear approximation of all possible 
scalar combinations of L x and L 2 . Note that as in the 
original BGK ansatz, we do not account for the possibil¬ 
ity of anisotropic scattering in turbulent fields. 

We now look for the power spectrum of a stationary 
solution of Eq. (12) using the approximation (13). The 
product (/ 1 / 2 ) can be expanded into eigenstates of J 2 
and J z and the power spectrum integral (1) corresponds 
to projections onto the state J = 0 with M = 0. The ro¬ 
tation term corresponding to the regular magnetic field in 
Eq. (12) does not contribute. Since the relaxation rates 
h'r(x) and v c (%) hi the approximation (13) are assumed 
to depend only on x = P 1 P 2 the relative rotation pro¬ 
portional to (L 2 + L 2 V 2 = L 2 is the only non-vanishing 
term. 


For identical momenta and rigidities, the relative dif¬ 
fusion should vanish completely. One way to see this is 
by considering a toy model where the initial configura¬ 
tion of CRs is homogeneous (Ahlers 2014). From this 
setup one can infer that the PSD at any time obeys 
d t f dp(/ 2 (t, r, p)) = 0. Therefore, in this case we must 
also have / dpzz r (l)L 2 (/ 2 (f, r, p)) = 0. Since this is in¬ 
dependent of the initial configuration we must therefore 
conclude that i/ r (l) = 0. 

A functional dependence can be obtained from a geo¬ 
metrical argument: two particles emitted locally under a 
small opening angle will separate with a reduced relative 
velocity Av ~ ^/2(1 — x) as long as both particles are 
experiencing a strongly correlated magnetic field. This 
leads to the ansatz v T (x) oc (1 — x) p with p = 1/2. To 
allow for a weaker or stronger ^-dependence, we have 
generalized this form to v r (x) oc (1 — x) p , with 0 < p < 1. 

Under these assumptions and approximations, the sta¬ 
tionary average Cg spectrum then obeys the set of equa¬ 
tions (see Appendix A) 

0 7, 1 1 f 

Qi$n = ^ ~2(Ck)k(k + 1)—-— / d xv T (x)Pg(x)P k (x ), 

k 

(14) 

with the effective dipole source term Qi = 
Kd, ndj'ri/ (6n) . Since the rhs of the previous equation 
is an expansion in terms of Legendre polynomials, we 
can easily derive the stationary spectrum as 


{Ci) 


3 Si j d x p d x ) 
2 £(£ + l) J V r (x) 


(15) 


With the assumed v r (x) oc (1 — x) p , this simplifies to 

= 3 fii ( p(p- 1) \ r(i -p)r(i + p-1) 

[t> 2 p \£(£+l) + J T(p)T(£-p + 3) • 

(16) 

In the limit i 1, Y(t + p — 1 )/T(£ — p + 3) — > £ 2p-4 , 
and so (Cg) oc £ 2p ~ 4 . For p = 1/2, this is close to the £~ 3 
behavior exhibited by the data (cf. Fig. 2). Figure 2 also 
shows examples of the power spectrum (15) for three 
models v r (x) oc (1 — x) p with p = 1/3, 1/2, and 2/3. 
While our ansatz v r (x) oc (1 — x) p is phenomenological in 
nature, a microscopic derivation of ^(x), in quasi-linear 
theory, will be presented in future work. 


5. CONCLUSIONS 

We have discussed the appearance of small-scale 
anisotropy from the propagation of CRs in local tur¬ 
bulent magnetic fields. We have shown that this effect 
can be understood as a phenomenon of relative diffusion, 
where CRs with similar arrival directions experience cor¬ 
related motion in the local magnetic environment. Our 
numerical simulation showed that the ensemble-averaged 
dipoles can be significantly larger than those expected 
from standard diffusion theory. In addition, the formal¬ 
ism predicts significant power in the ensemble-averaged 
small-scale multipoles. 

We have also provided an analytic framework that can 
describe the multipole power spectrum via an effective 
relative relaxation rate of multipole components in the 
evolution of the ensemble-averaged PSD products. This 
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formalism is a natural extension of the well-known BGK 
approach. Assuming that the relative relaxation rate has 
a simple power-law dependence on the relative opening 
angle of CR arrival directions, we arrive at a power-law 
spectrum very similar to that obtained from simulations 
and the observed CR data. 
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APPENDIX 

A. DERIVATION OF THE STATIONARY POWER SPECTRUM 

The decomposition of the PSD /(p) into spherical harmonics Y# m ( p) corresponds to a projection onto eigenfunctions 
| Em) of the angular momentum operator L. The products of PSDs that appear in the definition of the power spectrum 
in Eq. (1) are then eigenfunctions of \£imi£ 2 m 2 ) that can be decomposed into eigenfunctions | of the total 

angular momentum operator J = Li +L 2 . As noted earlier, the integral over arrival directions in Eq. (1) corresponds 
to a projection of / 1/2 onto states with J = 0 and M = 0 and £1 = £2 = £, which can be seen from the expansion into 
spherical harmonics 

\mt) = ^Ly, y L(pi)yUp2) = (-i) f ^^p,(p!p 2 ), (Ai) 

^ m 

following from the Clebsch-Gordan coefficients 

( _D V — m i 

(00|4toi£ 2 - m 2 ) = , • (A2) 

V2ti + I 

Now, the BGK ansatz (13) assumes that the scattering rates v r (x) and v c (x ) are only functions of x = P 1 P 2 and can 
thus be expanded into a sum over states with J = 0 and M = 0. The projection of Eq. (13) onto \00££) can thus only 
receive contributions from the projection of the product of PSDs onto J = 0 and M = 0 states, 

/r /2 = EE dpqQ'rs Ypq (pi )E rs (P2) ^ EE CLpqCL rs (00|p<3Ts)|00) = -E( 2 P+l)^(PiP 2 n- 

pq rs pq rs p 

Finally, the stationary solution (14) has the form Q\5n = J dxPe(x)K(x) and hence K(x) oc P\(x) = x. From this we 
immediately arrive at Eq. (15). 
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